o 

(N 
C 

a 

(N 
(N 



Model-based discrete relaxation process representation of band-limited 
power-law attenuation 



Sven Peter Nasholm 

Department of Informatics, 

(Dated: November 16, 2012) 



University of Oslo, P. O. Box 1080, NO-0316 Oslo, Norway 



Frequency-dependent acoustical loss due to a multitude of physical mechanisms is commonly mod- 
eled by multiple relaxations. For discrete relaxation distributions, such models correspond with 
causal wave equations of integer-order temporal derivatives. It has also been shown that certain 
continuous distributions may give causal wave equations with fractional-order temporal derivatives. 
This paper demonstrates analytically that if the wave-frequency uj satisfies SIl <C oj <C S1h, a con- 
tinuous relaxation distribution populating only SI £ [Ql , SIh] gives the same effective wave equation 
as for a fully populated distribution. This insight sparks the main contribution: the elaboration of 
a method to determine discrete relaxation parameters intended for mimicking a desired attenuation 
behavior for band-limited waves. In particular, power-law attenuation is discussed as motivated 
by its prevalence in complex media, e.g. biological tissue. A Mittag-Lefner function related distri- 
bution of relaxation mechanisms has previously been shown to be related to the fractional Zener 
wave equation of three power-law attenuation regimes. Because these regimes correspond to power- 
law regimes in the relaxation distribution, the idea is to sample the distribution's compressibility 
contributions evenly in logarithmic frequency while appropriately taking the stepsize into account. 
This work thence claims to provide a model-based approach to determination of discrete relaxation 
parameters intended to adequately model attenuation power-laws. 
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I. INTRODUCTION 

This paper concerns the determination of multiple re- 
laxation parameters. 

It is empirically observed that attenuation in biolog- 
ical tissue and other complex media such as polymers, 
rocks, and rubber often follows a power-law in frequency: 
afc(w) oc a;' ', with the exponent between and 2 ( Sz- 
abo and Wu (|2000f )). Such power-laws can be valid over 
many frequency decades. For acoustic modeling, time- 
fractional derivative wave equations have been shown to 
imply power-law at t enuat i on over wide frequency band s 
( Holm and Sinkusl (12010ft : iHolm and Nasholml (1201 il l). 
Such fractional wave equations can be obtained from 
linearized conservation of mass and momentum in com- 
bination with time-fractional constitutive relations con- 
necting stress and strain. Associate d nonlinear frac- 
tional wav e equations a r e pres ented in lPrieur and Holml 
(|201lh and|P ricur et all (|2012h while related linear wave- 
propagation i simulatiojis_^ire - j£parte^ 1 e i £ ii- ^^ 
and Ludwig dl995T): ILiebler et all (|2004l ): IWismerl ([20061 ); 
ICaputo et all (|201lh . 

More over, the multip l e rela xation mechanism frame- 
work of lNachman et all ( 1990T ) is widely considered ade- 
quate for acoustic wave modeling in lossy complex media 
like those encountered in medical ultrasound. It relies on 
thermodynamics and first principles of acoustical physics. 



The corresponding lossy wave equation for TV discrete 
relaxation mechanisms is a partial differential equation 
with its highest time derivative of order TV + 2. This 
model is in the following denoted the Nachman-Smith- 
Waag (NSW) model. 

Viscoelastic constitutive stress-strain models are gen- 
erally possible to convert into a Maxwell-Wiechert de- 
scription with springs and dashpots in parallel. The 
NSW model is linke d to fra ctional derivative modeling in 
iNasholm and Holml (|201lh . where a continuum of relax- 
ation mechanisms is assumed. The compressibility contri- 
butions were assumed to be distributed following a func- 
tion related to the Mittag-Lefner function. It was shown 
that the wave equation corresponding to this distribution 
is identical to the fractional Zener wave equation. Ac- 
tually, a Maxwell-Wiechert description of the fractional 
Zener str ess-strain relation was thereby implicitly ver- 
ified. In lAdolfsson et al. I (|2005l ). a very large number 
of weighted Maxwell elements evenly distributed in the 
linear frequency scale are shown to give the same stress 
response as a fractional order viscoelastic model. We 
also note that the rheological fractional spring-pot ele- 
ment was i nterpreted in te r ms of weighted springs and 
dashpots in lPapoulia et all ( 20101 ). 



Kelly and McGough (|2009h demonstrated that hierar- 
chical fractal ladder networks of springs and dashpots can 
lead to power-law attenuation in a low-frequency regime. 
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This approach however requires a large number of degrees 
of freedom which makes parameter fits cumbersome. 

Band- limited fits to power-law acoustic attenuation for 
relaxation models with N = 2 and 3 are exemplifie d by 
iTabei et all (120031 ) and lYang and Cleveland! (j2005T ). In 
the latter, one of the mechanisms is assumed to be of very 
short relaxation frequency, thus representing a thermovis- 
cous component. The two other mechanism relaxation 
frequencies, the corresponding two compressibility con- 
tributions, and the compressibility of the thermoviscous 
component are determined through numerical minimiza- 
tion of the difference between the resulting attenuation 
and the desired power-law. 

For a large number of modeled relaxation mechanisms, 
such numerical optimization of the parameter fit turns 
very intricate. 

Recently published works which stress the need for 
straightforward and accurate determination of discrete 
relaxation repr esentations that gene rate power-law atten- 



uation, are e. 



Roitner et al. | (|2012t ) (see Section II) and 



iTreebv et all (j2012T) fsee Sectio n II. C). See also the intro 



duction in lLiebler et ~d\ (120041) . 

The present paper establishes a systematic model- 
based method to choose the compressibility contribution 
and the relaxation frequency for each of N relaxation 
processes, to model power-law attenuation over a given 
wave frequency band. This discrete parameter selection 
method is based on the continuous distribution of relax- 
ation processes previously shown to result in fractional 
Zener wave equation, which in turn generates 3 distinc- 
tive p ower-law attenuation regimes ( Nasholm and Holml 
(2011)). The present work also studies analytically the 
effect of letting this distribution cover only a limited fre- 
quency band. 

This paper is organized as follows. The Theory Sec- 
tion [IT] first reviews and partially extends relevant NSW 
theory, fractional Zener model considerations, and the 
link between those. Then it presents new developments 
related to band-limited continuous and discrete relax- 
ation distributions, as well as the related attenuation for 
waves of frequencies within and outside the populated 
relaxation bandwidth. Section III. El is central because 
it describes how to select discrete relaxation process pa- 
rameters to get approximate power-law attenuation over 
a given wave-frequency band. Section IIIII provides two 
numerical examples exemplifying how discrete relaxation 
parameters may be determined in order to attain power- 
law attenuation. Discussions and conclusions are given 
in Section [TV] 



II. THEORY 

A. Relaxation processes modeling within the NSW 
framework 

1. Discrete relaxation distribution 

The NSW model of multiple discrete relaxation pro- 
cesses results in the frequency-domain generalized com- 



pressibility ( Nachman et all ( 19901 )) 
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where the mechanisms v = 1 . . . N, have the relaxation 
times T\, . . . , tjv and the compressibility contributions 
K\, . . . , Kjy. In the present work il„ = 1/t v denotes the 
process relaxation frequency. 

The frequency-domain generalized compressibility is in 
the following simply denoted compressibility. In some 
branches of science it is instead called complex compli- 
ance J*(oj), being defined as the ratio between strain and 
stress: k(lu) = e(w) / 'cr(w). It is hence directly related to 
the constitutive stress-strain relation. 



2. Continuous relaxation distribution 



Following iNasholm and Holml ([20111 ). a representa- 
tion of Eq. ([1} when considering a continuum of re- 
laxation mechanisms distributed in the frequency band 
fl £ [flh , ^h] with the compressibility contributions de- 
scribed by the distribution ft„(il) becomes 

kn(w) = Kq — iuj [ — - dSl. (2) 

Letting the integral go from JIl = to fin = oo, and 
instead incorporating any possible relaxation distribution 
bandwidth limitation of k v ($X) into the distribution, we 
define = H(Q - n L )H{Q H - il)K u (Si), where H(Q) 

denotes the Heaviside step function. Then the integral 
@ is a Stieltjcs transform. By application of the Laplace 
transform relation, which for real lo, real t, and real f2 is 
valid for SI > 0: 



/>oo 

£ t {#(i)e-^}(il) 4 / H(t)e 
Jo 



iujt „ — fit 



dt = 



1 



(3) 



from the t domain to the Q domain, the compressibility 
@ may be written as 

kn(w) = K -iuj n„(n)^J^ ff(i)e" nt e- iwt dtjdn, 

(4) 

which by virtue of Fubini's theorem may be written as 



kn(w)= kq — ioj 



H(t) / K v (Q)e- ut dCl}dt 



kq - iuF t {H{t)£n {K„(n)Xt)}(w), 



(5) 
(6) 



where the Fourier transform from the t domain to the cu 

/oo 
e~ iult f{t) dt. 
-OO 
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Another equivalent representation of the compressibil- 
ity © is 

kn(w) = kq — w / — — dS2 — zu; / — — di 2, 



model afe(w) may hence be constructed using the inverse 
transform recipe 



o n 2 + ui 2 



o n 2 + 



(7) 



which is found by multiplying both the numerator and 
the denominator of Eq. by (il — ius). 



3. Attenuation and phase velocity 

The conventional decomposition of the frequency- 
dependent wavenumber k(u>) into its real and imaginary 
parts, gives the phase velocity c p {w) = {k(uj)} and 
the attenuation a^iuS) = — Q{/c(u;)}. 

Combining the definition of k(lu) = e(uj)/o~(u)) with 
the linearized conservation of mass and momentum (see 
iNasholm and Hohnl ( 201 lh and the references therein for 
details) gives 



k 2 (oj) = oj 2 pqk(oj). 



(8) 



In general, the attenuation and the phase velocity are 
thus given from the dispersion relation above as 



c» = w/B {k} = P - 1/2 /5R { v^R} ■ 



and 



(9) 



Analyzing this expression further, we see that in the 
small-attenuation regime where the kq part of (J7]) domi- 
nates over the imaginary part infers 



\AnM 



= \ K -U z 



n 2 + w 2 



dO - 



n 2 + uj 2 



do - 



n 2 



(10) 



Then the attenuation ctk is given from the approximation 



AS 



D 2 



(11) 



= Au 2 / / ft(Q)e- nt cos(wt) dt dn 
Jo Jo 

= A^ c {£ n {«*(n)}(*)}M> 



(12) 



where T c denotes the Fourier cosine transform from the 
t domain to the u domain, and £q denotes the Laplace 
transform from the f2 domain to the t domain. The in- 
troduced A is a frequency-independent scalar. The inte- 
gral o f Eq. (jlip is a Widder potential transform ( Widderl 
(11966m. 

Using Eq. © and (|12l) under the small-attenuation 
constraint, the relaxation process spectrum k v (Q) cor- 
responding to some frequency-dependent attenuation 



«„(fi) = A ■ £7 < T l 



Qfe(w) 



(*) (n). 



(13) 



Expression (| 13[) is si milar to what was fou nd using an 
approach reported in IVilenskv et all ( 2012[ ). however a 
formula eq uivalent t o (1111) was used already in e. 0. Paulv 
and Schwan (|l97ll ). The small-attenuation assumption is 
probably reasonable for compressional wave propagation 
in biological tissue. By contrast, for shear-wave propaga- 
tion, the attenuation i s generally much more pronounced 
()Szabo and Wul (|2000h 1. 

For a discrete set of relaxation mechanisms under the 
small-approximation assumption, we see from Eq. (|lip 
that the total attenuation is just the sum of the contribu- 
tion from each mechanism: afe(w) = J2u=i a ^( w )j where 
the contribution of mechanism v is given by 



,( W ) = 



n?,' 



(14) 



Such behavior is especially well documented in air fBass 
e±_aL (119951)")" and in sea water (lAinslie and McColm! 
(|1H98T)1. 



B. The fractional Zener wave equation 

As a consequence of t he fractional Zener model stress- 
strain relation (see e.g. iBaelev and Torvikl (1983)), the 
frequency-domain fractional Zener compressibility is ob- 
tained from the ratio e(uj)/o~(uj): 



A 1 + {r^f 

KZ[UJ) = K —— — -. 

1 + (r ff iu) a 



(15) 



Due to thermodynamic constraints, < a < 1 and f3 < a 
( Glockle and Nonnenmacherl ( 199lh l . however the case 
a = j3 is the most well-behaved from a physical point 
of view iRossikhin and Shitikoval ( 200lh . Insertion of the 
compressibility (fl"5"j) into the dispersion relation @, gen- 
erates the fractional Ze ner dispersion relation f Holm and 
NasholmHoilJ)) 



k 



UJ 2 1 + {Tjjjf 

cl 1 + {T a iu}) a 



(16) 



This is the spatio-temporal frequency domain representa- 
tion of the five-parameter fractional Zener wave equation 
which by inverse transform hence becomes 



vV 



1 d 2 u 
dUt 2 



dt a 



V 2 u 



t[ dP+ 2 u 
^W+ 2 



= 0. 



(17) 



The a = (3 variant of the fractional Zener compressibil- 
ity (|15|) combined with Eq. (HQ) results in three di stinct 
attenuation power-laws ( Holm and Nasholml ( 201lh ): 



io l+a low-frequency regime, 
otk oc ^ u) x ~ a l 2 intermediate frequency regime, (18) 



a,' 



1-a 



high-frequency regime. 
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C. Connecting the NSW and the fractional Zener models 

1. The continuous relaxation distribution Ac£,(f2) 

Provided that the linearized conservations of mass and 
momentum are valid, and given that the NSW com- 
pressibility «n(w) of Eq. (|6]) is equal to the fractional 
Zener compressibility kz(w) of Eq. (|15p . the dispersion 
relations from Eq. ([5]) are also equal. Because the dis- 
persion relation is a spatio-temporal Fourier representa- 
tion of the wave equation, kn(w) = kz(cj) thus implies 
that the NSW wave equation becomes equal to the frac- 
ti onal Zener wave equation (|17[) . As brought forward 
in lNasholm and Holml ( 201 lh . there exists a continuous 
distribution k' u (CI) of NSW relaxation processes which 
implies kn(^) = ^z(w) when a = /3: 



<(n) = 



1 k (t£ - Tf)fi a_1 sin(cnr) 
7T (r a n) 2a + 2{T a Q) a cos(aTr) + 1 ' 



(19) 



This distribution is the inverse Lap lace transform of a 
Mittag-Leffler related expression, see lNasholm and Holml 
(1201 lh :or details. Note that this link between the frac- 
tional Zener and the NSW models is valid also outside 
the small-attenuation regime 9{/e} <C 5ft {k}. 



2. k^(O) power-law regimes 

Below follows an analysis of the relaxation time dis- 
tribution k' v {Q) of Eq. (fH)j) as it evolves depending on 



1 k (t* - r^n"- 1 sin(aTr) 
7T (T a V,) 2a + 2{T a Vl) a cos(a7r) + 1 
C L ■ n a ~\ for ftr a < 1 

a • rr 1 , for n Ta « i 

C* H ■ n- a ~\ for 1 < n Tr7 , 



where the frequency-independent constants are: 
a k (t^ ~ rf)sin(a7r) 



(20) 



Ci 



a. «o(t° - T e ")sin(Q!7r) 

27rr°(l + cos(a7r)) 
a K (r° - r e a )sin(o!7r) 

__2a 



and 



(21) 



Due to the restriction < a < 1, we hence note that 
depending on the value of the product T a Q, the distribu- 
tion K' v (fl) may have the form of a power-law k! v oc ti d , 
where -2 < d < 0. 



i. Formal analysis of band-limited power-law relaxation 
distribution 

The case of a relaxation process continuum populat- 
ing the region ft <G [Hlj^h], for example where K„(f2) 
is adequately be approximated by a power-law (see Sec- 
tion HTC2] above) 



C 




ft" 



fte [n L ,n H ] 

otherwise. 



(22) 



where —2 < d < 0, is now investigated. Distributions 
of such form are related to broken power-laws, often de- 
noted Pareto distribut ions, and are co mmonly observed 
in nature and society ( Newman! ( 20051 )1. 

For Kv(fl) = C ■ ft d , the generalized compressibility © 
may be written 



kn(w) = kq — iuiC 



ft" 



ft 



■ dft 



ft 



dft 



(23) 



which in view of Eq. (J7]) and lBateman and Erdelvil (|1954 ) 

[Eq. (21)], becomes 



«n(w) = K + 



ft 



d-l 
H 

- d 



l-d 



iCuj 



(24) 



where $(a;, ^ , w) is the analytic continuation of the Lerch 
transcendent. In Section Til. D. 2 1 below, the compressibil- 
ity (f2"4")l is further analyzed. First a continuous band- 
limited power-law relaxation spectrum is considered and 
the corresponding generalized compressibility is deduced. 
Then three different wave-frequency regimes are consid- 
ered for any band- limited relaxation relaxation K„(ft) 
and the corresponding attenuation functions are brought 
forth. For v = 1, the function §{x,v, u) is related to 
the analytic continuation of the Gauss h ypergeometric 
function as: (|Bateman and Erdelvil (jl954T ). Eq. (1.10)) 



l 2 F 1 (l,u; 1 + u;x), 



(25) 



to which there is a dire c t conn ection to fractional calculus 
because ( Samko et all ( 19931 )) 



iFi(a, b;c;x) 



T(a)a; 1 



„b-l 



T(b) 



(1 



(26) 



D. A band-limited continuum of NSW relaxation processes 

In the following, the case of a continuum of relax- 
ation processes populating a limited frequency band ft G 
[ftL j ftfl] is further explored with respect to the for three 
different wave-frequency regimes. 



For practical purposes when u> is far from the relax- 
ation process cut-off frequencies fti, and fin, it is how- 
ever not necessary to evaluate the Lerch transcendents 
of (|24[) to find k(ui) and hence determine the frequency- 
dependent attenuation and phase velocity. The three rel- 
evant frequency regimes are considered in the following 
two subsections. 
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2. The Ql -C ui -C frequency range for any k„(w) 

In case the wave-frequency is within the frequency 
band of the relaxation processes, in the integrals of ([7]) 
we change the integration limits from into 51l and from 
oo into fin and in addition make the variable change 
u = fl/ui to result in 

kk(u) = K - U) - Au-iui du. 

I S!l u z + 1 JIIl + 1 

(27) 

For the analyzed case ftL *C w ftij, the integration 
limits then become Ol/w - ^ and ftn/w — > oo. Then 
performing the variable change back into 17 = uui makes 
it clear that under the given conditions the effective com- 
pressibility thus becomes equal to the compressibility re- 
lated to the fully-populated relaxation frequency band 
([7]). A similar analysis may also be done directly in ©. 
The fractional Zener wave equation (|17[) and its asso- 
ciated attenuation and ph ase velocity expressions (see 
iNasholm and Holml (l201ll )) are consequently valid also 
for the band-limited distribution of relaxation processes, 
as long as the wave frequency ui is much lower than S1h 
and much greater than fti,. 

Interestingly, the analysis above prompts the conclu- 
sion that for a band- limited (Ol <C oj <C J7h) power- 
law relaxation distribution K„(ft) oc ft d , which results in 
the compressibility «n(w) given by (|24p . there is a di- 
rect relation to the attenuation power-law valid within 
ft L < w < n H -. 

First, consider the low-frequency power-law regime of 
the distribution fi^(ft) in Eq. (|19p being populated in the 
limited relaxation frequency band. Then k'„(^) k ^ d 
which gives d = (a — 1) e [—1,0], because a <E [0,1]. 
By virtue of Eq. (|18|) this corresponds to the attenuation 
power- law ak(t-o) oc u/' with the exponent r/ £ [1 + a] <^4> 
77 £ [1,2]. Second, consider instead the high-frequency 
regime of the distribution n' v (Q,) being populated in the 
limited relaxation frequency band. This leads to ftj,(ft) oc 
fl d with d = (-a - 1) £ [-2, -1]. 

Hence a band ftL "C cj <C ftfj populated with relax- 
ation frequency contributions ft„(ft) oc ft d gives a band- 
limited power-law attenuation a.k{ui) oc w d+2 both for 
d £ [-2,-1] and for d £ [-1,0]. 

3. The ui -C fi/. < 51h and 51/. < fin -C ^ frequency regimes for 

any k, v (ui) 

As pointed out in lNasholm and Holml ()201ll ). the case 
of the wave-frequency being much lower than the popu- 
lated relaxation process frequency band ft € [fti,, ffipK an 
analysis similar to the LF limit of iNachman et all (|l990f ) 
results in ctk oc us 2 . By the same token, for the wave- 
frequency being much higher than the populated relax- 
ation proportional to frequency band, the attenuation is 
frequency- independent . 

Within these frequency regimes, the effective wave 
equation is hence the same as for a single discrete relax- 
ation mechanism with the compressibility contribution 



being proportional to the impulse <5 (O — (ftL + ftH)/2). 
This is similar as for a medium with one single discrete 
NSW relaxation process. 

E. N discrete NSW relaxation processes to get 
band-limited attenuation power-law 

Based on the developments in Section III.D1 here a 
method is elaborated to determine the ft„ and n u pa- 
rameters of a discrete set of TV relaxation processes in- 
tended to result in an appropriate attenuation power-law 
otk(oj) oc u/ 1 for waves within the bandwidth ui £ [ftL, ftii]- 
The power-law exponent satisfies rj € [0, 2]. 

First, based on (|18p . the r\ exponent determines 
whether a low, high, or intermediate frequency model 
is applied. Based on this, t g is then set so that \/r a 
is either much higher or much lower than the frequency 
region of interest [S1l,^h]- Subsequently, the relaxation 
process frequencies ftj, are sampled within [ftL, ^h], equi- 
spaced in the log ft domain hence giving 

N-u v-\ 

n v = ft™- 1 ft™- 1 . (28) 

Thereafter, the compressibility contribution of 
each process is decided from n' v (£l v ) of (TT9|) . Alterna- 
tively, the relevant power-law approximate expression of 
(|2U)) may be chosen. Finally, the T e parameter is adjusted 
to achieve the attenuation oik (w re f ) = a re f at some appro- 
priate reference frequency us lc f € [ftL, flu]. 

For calculation of the resulting compressibility kn(u), 
the discretized approximation of the integral ([5]) must 
take into account the uneven stepsize 

Aft, = ft, (l - (ft L /ft H ) 1/(Ar_1) ) . (29) 

The kn(w) estimate kn(w) which approximates ([2} for 
the chosen discrete relaxation process parameters is 
hence 

kn(w) = kq — iui — ^- — — Aft,,. (30) 
^— ' ft„ + IUS 

v—l 

From this compressibility, the attenuation and phase ve- 
locity may be calculated in the conventional manner us- 
ing Eq. The attenuation resulting from this relax- 
ation process parameter decision approach is explored in 
the following 2 numerical examples. 

III. NUMERICAL EXAMPLES 

A. Power-law a k oc w 1 ' 1 for / £ [100 kHz, 30 MHz] 

Below follows an explicit example of modeling an at- 
tenuation power-law ak oc ui between 100 kHz and 
30 MHz using discrete relaxation processes. Such atten- 
uation is relevant for medical ultrasound imaging. The 
chosen medium properties, which are l isted i n Table HI arc 
the same as in lYang and Cleveland! (l2005h . The atten- 
uation resulting from the parameter selection approach 
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TABLE I. Medium parameters fo r the attenuation power-law 
fit of Section IITOI similar to the I Yang and Cleveland! (|2005h 
parameters. 



Equilibrium speed of sound, Co (m/s) 
Density, p (kg/m 3 ) 

Zero-freq. compressibility, /-co = -y— (Pa -1 ) 4.0158 

c qP0 

Wanted attenuation at 1 MHz, Qo (dB/MHz/cm) 
Wanted attenuation power-law exponent, r\ 
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FIG. 1. (Color online) Attenuation as a function of frequency 
for TV o-relaxation mechanisms determined by the method 
proposed in Section III. El aiming at constructing a power-law 
valid within 100 kHz and 30 MHz. Also displayed are the 
lYang and Cleveland! (|2005h N = 3 attenuation and an ideal 
power-law oc u 1 ' 1 (thick solid line). Table [T] lists the medium 
properties. Top pane: Attained attenuations. Bottom pane: 
Relative difference between the attained relaxation attenua- 
tion functions and the power-law. The horizontal axes rep- 
resent wave-frequency. For visualization convenience, each 
attenuation is normalized in order to make the minimum and 
maximum relative differences equal in magnitude within the 
relevant frequency interval. The Yang and Cleveland attenu- 
ation is included both with normalization (thick dashed line) 
and without (thin dashed line, only included in the bottom 
pane). The attenuation is shown for TV = 1 (thin dashed 
line), TV = 2 (thick dash-dotted line), and TV = 3 (thin solid 
line). The relaxation parameters fl„ and k v are displayed in 
Table IH 



proposed in Section III. El was applied for TV = 1 , . . . , 4 
relaxation mechanisms using S7l = 100 kHz and Oh = 
10 MHz. The consequent attenuation functions as calcu- 
lated from the NSW compressibility in Eq. ([T]) are dis- 
played in Fig. [TJ where also a pure oek{ui) oc ui law is 
pl otted, as well as the Yang and Cleveland attenuation. 
In I Yang and Cleveland! (I2005IV a numerical least squares 



TABLE II. Relaxation process frequencies S7„ and compress- 
ibility contributions k v corresponding to the attenuation func- 
tions displayed in Fig[Tj The parameters were attained using 
the Section III. El method for TV — 1,2,3 relaxation mecha- 
nisms. 
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scheme was used to fit two relaxation terms and a thermo- 
viscous component to an approximate oc of 1 power-law 
using Eq. (|14[) . The relaxation frequency £l v was pre-set 
for the thermoviscous component, therefore leaving two 
17„ and three k v parameters to be determined by the nu- 
merical scheme. The discrete relaxation frequencies and 
compressibility contributions arc displayed in Table |TTJ 

A fit for the same attenuation law to the fractional 
Zener model and to the correspon ding continuous com- 

ressi bility function is provided in iNasholm and Holm! 

201 ll ) where Table II lists the parameters a, r e , r a , and 

K . 



B. Power-law ct k oc w 1 - 1 for / G [100 kHz, 1 GHz] 



In the following, the approach suggested in Section llI.EI 
is applied in a similar way as in Section rill.Al to construct 
power-law attenuation within the wider wave-frequency 
band / <E [100 kHz, 1 GHz]. This frequency band cov- 
ers more or less all frequencies of conventional pulse-echo 
medical ultrasound imaging and microscopy. The prop- 
erties of discrete relaxations with VL V evenly distributed 
in log n were determined for distribution sets with 1 to 
7 mechanisms. The attenuation functions corresponding 
to each relaxation mechanism set are displayed in Fig. [5J 



IV. DISCUSSION AND CONCLUDING REMARKS 



As laid out in Section III.D1 the present paper demon- 
strates that for a continuous distribution of relaxation 
processes with the compressibility contribution distribu- 
tion Kv(Q) populating the relaxation frequency range 

G [Qlj^h]) the effective wave equation for the wave- 
frequency oj satisfying 17 l <C u> *C 17h is the same as for 
the relaxation processes populating the whole 17 G [0, oo]. 
This work thus supports the intuitive conjecture that for 
a distribution of relaxation processes covering all frequen- 
cies, it is the 17^ within the wave- frequency bandwidth 
that mainly contribute to the attenuation. The wave may 
thus be seen as probing the medium around the relax- 
ation frequencies within and close to the wave-frequency 
bandwidth. 

When instead cj <C S7l < 17h or S7l < 17h £*;, the ef- 
fective wave equation becomes the same as for a single dis- 
crete NSW relaxation mechanism at 17„ = (17l + 17h)/2. 
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FIG. 2. (Color online) Attenuation as a function of frequency 
for TV relaxation mechanisms determined by the method pro- 
posed in Section III. El aiming at constructing a power-law 
oc oj ' (thick solid line) valid within the wide interval between 
100 kHz and 1 GHz. The horizontal axes represent wave- 
frequency. Top pane: the attenuation functions for N — 1 
(thin dashed line), N — 2 (thick dashed line), and N = 4 
(thin solid line). Middle pane: Relative difference between 
the attained relaxation attenuation functions and the power- 
law, with each attenuation normalized in order to make the 
minimum and maximum relative differences equal in magni- 
tude within 100 kHz and 1 GHz for N = 4 (thin solid line), 
N = 5 (thin dashed line), TV = 6 (thick dashed line), and 
N — 7 (thick dash-dotted line). Bottom pane: Same as the 
middle pane, however with the normalization taking the more 
narrow frequency interval between 500 kHz and 0.2 GHz into 
account. 



All attenuation models that may be written as a super- 
position of NSW relaxations are causal. 

In particular this work considers the Mittag-Lcfflcr 
function related distribution case k' v (U) as given in ([2T))) . 
For a fully populated range Q G [0, oo] this has previously 
been shown to produce the four-parameter fractional 
Zener wave equation prompting three distinct frequency 
power-law attenuation regimes. If this relaxation mecha- 
nism distribution instead populates only Q G , ^h] , as 
shown in Section Hi. D. 21 the 4-parameter fractional Zener 
wave equation is still valid as long as the wave-frequency 
u) satisfies SIl <C oj -C ^h- Actually, the distribution 
n' u {Vl) may, as shown in Section Til. C. 21 for three differ- 
ent f2 regions be appropriately approximated by oc oj d , 
with — 2 < d < 0. These regions are linked to the atten- 
uation power-law frequency regimes of the 4-parameter 
fractional Zener wave equation. 

The insights of Sections III. C I and III. D I provide the basis 
upon which a proposed method for selection of discrete 
relaxation distribution properties relies. The discrete re- 



laxation parameter parameter selection approach is laid 
out in Section III. El which is intended to be a key ad- 
vance of the current paper. One may hence avoid ap- 
plying numerical minimization methods to determine the 
relaxation parameters, which becomes cumbersome espe- 
cially when there is a large number of mechanisms. This 
proposed selection of discrete parameters is model-based 
primarily in the sense that for a chosen process relaxation 
frequency £L Vl the compressibility n v is decided from the 
distribution k' v {Q, u ) given in (|20|) . The 27V-parameter 
determination task to decide the fl v and k v values for 
a model of N relaxation mechanisms thus becomes an 
A^-parameter determination problem. However for deter- 
mination of the N relaxation frequency parameters, the 
straightforward choice to distribute £l u logarithmically 
equi-spaced within which il„ G [f^L, Oh] is proposed here. 
The limits Ol and f2u are set similar to the maximum and 
minimum wave-frequencies where the attenuation model 
is to be valid. 

As power-law attenuation is commonly encountered for 
large wave-frequency intervals in complex media, this 
work puts emphasis on the determination of discrete 
relaxation parameters to yield power-law attenuation 
within a given frequency interval. The theoretical con- 
siderations of Section [II] are supplemented by two numer- 
ical examples in Section [TTTl For the first example, where 
the attenuation is displayed in Fig. [TJ we note that for 3 
relaxation the chosen weighted mechanisms give rise to 
an attenuation differs by less than 14% relative to the 
pure power-law a.k{w) oc uj 1a within the whole frequency 
region of interest / G [100 kHz, 30 MHz]. For two mech- 
anisms, t he highest relative dif f erence is 63%. By con- 
trast, the lYang and Cleveland! 1)20051) weighted sum of 
three relaxations corresponds to a maximum relative dif- 
ference in a.k{oj) of 68%. When normalizing the Yang 
and Cleveland attenuation, the maximum difference is 
reduced to 29%, which is however still more than dou- 
ble the 14% limit attained when setting the parameters 
following the prescription suggested in the current paper 
for the same number of relaxation mechanisms. Judging 
from the maximum relative absolute difference between 
the resulting ctkito) and the pure power-law with the ex- 
ponent 1.1, the proposed method hence is more advanta- 
geous. 

Inspection of the resulting attenuation curves of Fig. Q] 
hints the possibility of attaining less relative difference 
between ak(u)) and the power-law by straightforward ad- 
justment of the width of the populated relaxation fre- 
quency region. This possibility is now briefly explored. 
For low N, it seems appropriate to make the populated 
region smaller, while for larger A^ the relative difference 
between attained afe(w) and the wanted power-law in- 
creases close to f^L and fin hence suggesting to make 
the populated frequency region wider. For N = 2, 
shrinking the populated relaxation frequency range from 
n G [0.1, 30] MHz via [0.2, 15] MHz down to [0.3, 11] MHz 
decreases the maximum relative difference from 63% via 
39% down to 28% (see Fig. [3]), which is actually about 
the same as was attained for the normalized Yang and 
Cleveland fit for three relaxation processes. 

The second numerical example for which the resulting 
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FIG. 3. (Color online) Relative difference between the at- 
tained attenuation and the wanted power-law oc u for 
N — 2 mechanisms, as attained when the relaxation fre- 
quencies are fii, 2 = {0.1, 30} MHz (solid line), {0.2, 15} MHz 
(dash-dotted line), and {0.3, 11} MHz (dotted line). All these 
attenuation functions are normalized so that absolute max- 
imum and absolute min imum of the relative diffe rence are 
equal. For reference, the I Yang and Cleveland! (|2005f ) attenua- 
tion is also included, both in the original form (thick dashed 
line) and in the adjusted normalized form (thin dashed line). 
The vertical lines indicate frequencies of the explored relax- 
ation mechanisms. 



attenuation functions are displayed in Fig. [3J illustrates 
that that for as few as five discrete relaxation processes, 
an adequate power-law differing by less than 11% from 
the wanted power-law may be constructed within the very 
wide frequency band / € [100 kHz, 1 GHz]. 

The fractional Zcncr compressibility of Eq. (|T5j) is for 
a = f3 similar to the Cole-Cole expression for complex 
dielectric permittivity, which is empirically shown to be 
valid in a variety of complex media. A single discrete 
acoustical NSW relaxation mechanism in the dielectric 
permittivity context corresponds to a single Debye term. 
The adjustment and weighting of multiple discrete De- 
bye terms to emulate Cole-Cole dielectrical beh avior is 
treated e.g. in lRekanos and Papadopoulosl (l2010h. where 
a Pade approximant approach is applied, in lTofighil ( 2009f ) 
where an error-minimization method is employed choose 
the relaxation times at whi ch the Cole-Cole rel axation 
distribution is sampled, in iKellev et al ] (120071 ) where 



the relaxation frequencies are found using a nonlinear 
method and the weights are found using a linea r least 
squares approach, as well as in lClegg and Robinson! (|2010h 
where a genetic algorithm is applied to find the multi- 
ple Debye parameters. The parameter determination ap- 
proach brought forward in the present paper should also 
be tested out for selection of a discrete set of Debye terms 
to emulate the Cole-Cole dielectric permittivity. 

The model applied in iBerkhoff et al. I (119961 ) was con- 
tcxtualizcd in the discussion of Nsholm and Holm (2011). 
It was found that it corresponds to band-limited continu- 
ous relaxations with compressibility contributions given 
by K^ifl) cx fi . 

W e note that based on experimental evidence. Jongen 
et al. (| 19861 ) suggests the attenuation for example in beef 
liver to be cx oj 2 below a certain cut-off frequency uj c 
and oc w at higher frequencies. Such behavior is attained 
by the band-limited continuous relaxation distribution 



framework of the present work when the lowest frequency 
f^L of the populated relaxation frequency band is equal 

tO LU C . 

The N relaxation parameter determination method of 
Section III. El corresponds to mapping a time-fractional 
wave equation into an integer-order one of highest order 
N + 2 with the intent of the mapping being adequately 
valid within a given wave-frequency band. Because 
fractional-order differential wave equations require extra 
care in simulations, the method suggested here hence 
has a potential to facilitate numerical wave-propagation 
calculations. Comparison to previously suggest e d such 
numerical schemes as e.g. IWismer and Ludwid (Il995l ): 
ICaputo et al\ ( 201 lh . is an appropriate future work con- 
nected to this paper. 

From a practical point of view the method of this paper 
to determine discrete relaxation parameters is appealing 
because the anomalous physics of fractional-order differ- 
ential equations is converted to a differential equation 
with an finite set of higher-order integer derivatives. This 
conversion is tailored for the frequency bandwidth of in- 
terest. For such conversions to be valid for all frequencies, 
an infinite number of terms are needed in the integer- 
order differential equation. Thence the narrowing of the 
frequency region is traded off into the convenience of hav- 
ing a finite highest derivative order. 

Finally this work calls upon further experimental ver- 
ifications and more profound theoretical connection be- 
tween the relaxation models and fundamental physical 
properties of complex materials. Hopefully the links be- 
tween fractional calculus, observed attenuation behav- 
iors, the stunningly common power-law patterns of na- 
ture, and the micromechanical structure of matter are to 
be further clarified. 
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